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—1. . The divergence of the heat conductivity in the thermodynamic limit is in- 

cn . vestigated in 2d-lattice models of anharmonic solids with nearest-neighbour 

^^ \ interaction from single-well potentials. Two different numerical approaches 






based on nonequilibrium and equilibrium simulations provide consistent indi- 
cations in favour of a logarithmic divergence in "ergodic", i.e. highly chaotic, 
dynamical regimes. Analytical estimates obtained in the framework of linear- 

^L , response theory confirm this finding, while tracing back the physical origin of 

'^ \ this anomalous transport to the slow diffusion of the energy of long- wavelength 

O \ effective Fourier modes. Finally, numerical evidence of superanomaious trans- 

^ ' port is given in the weakly chaotic regime, typically found below some energy 

CJ ■ density threshold. 
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I. INTRODUCTION 

The study of heat conduction in models of insulating solids is a long-standing and de- 
bated problem. In 1929 R. Peierls provided the first convincing theoretical explanation of 
this phenomenon, relying upon the hypothesis that lattice vibrations responsible for heat 
transport in insulating solids can be described as a diluted gas of quantum quasi-particles, 
the phonons [|l|. Peierls' model extends Boltzmann's kinetic theory of the classical ideal 
gas to the phonon gas, by substituting collisions in real space with Umklapp processes in 
momentum space. Due to the quantum character of this approach, Umklapp processes 
were theoretically justified as a perturbative effect, epitomizing the scattering of phonons 
originated by intrinsic features of real solids: disorder and nonlinearity. The main achieve- 
ment of Peierls' theory is the prediction of the proportionality relation between the thermal 
conductivity k, and the specific heat at constant volume cy: 

k{T) = yv{T)vX (1) 

where v and A are the average velocity and the mean free path of phonons, respectively. 
Using Debye's formula, the dependence of k on the temperature T was found to agree with 
experimental observations in the low-temperature regime. Nonetheless, at least in the high 
temperature regime, there is no reason why a purely classical model of an insulating solid 
should not reproduce the correct behaviour of /t(T) and Fourier's law 

J{x) = kVT{x) (2) 

where J{x) and T{x) are the heat flux and the temperature field, respectively. To our 
knowledge, the first attempt of its theoretical justification on the basis of a microscopic model 
of a solid was made in a seminal paper by Rieder, Lebowitz and Lieb . They considered 
a chain of A^ coupled harmonic oscillators in contact at its extrema with stochastic thermal 
baths at different temperatures. They prove that in this simple model of a purely harmonic 
solid thermal conductivity diverges in the thermodynamoc limit as k,{N) ~ A^. Consistently, 
the temperature profile in the bulk of the chain is fiat. The physical interpretation of these 
results is that harmonic lattice waves, i.e. Fourier's modes, freely propagate through the 
lattice, thus contributing to a ballistic rather than diffusive transmission of heat. The 
same result was obtained numerically for another integrable model, the Id Toda lattice 0. 
Its normal modes, the so-called Toda solitons, are nonlinear waves localized in space and 
also freely propagating through the lattice. More generally, any integrable Hamiltonian 
system is expected to be characterized by anomalous transport properties since its normal 
modes, in the perspective of kinetic theory, are equivalent to a "gas" of non-interacting free 
particles. Following Peierls, the addition of generic disorder and nonlinearity should overtake 
these anomalies by introducing effective diffusive behaviour. Numerical [Q and analytical 
estimates performed for the mass-disordered harmonic chain still yield anomalous transport, 
k{N) ~ A^^/^, although the temperature gradient was found to be nonzero ||^. Accordingly, 
disorder is not sufficient to guarantee normal transport properties. It remains nonlinearity, 
that since the 70 's has been widely studied also as a source of chaotic behaviour. In fact, 
despite its deterministic nature, chaos was shown to yield good statistical properties, namely 
ergodic behaviour, for most of the thermodinamically relevant observables. Accordingly, 



many attempts were made with models of chains of nonhnearly coupled oscillators, that did 
not achieve convincing conclusions (e.g., see [0-^]). Finally, Casati et al. [jlO| proposed the 



ding-a-ling model that exhibits normal heat conductivity. It is worth mentioning that the 



same result has been confirmed one decade later by Prosen and Robnik [^ for the ding-dong 
model - a variant of the ding-a-ling one. The great merit of the former paper has been to 
point out the possibility of obtaining normal conductivity in a nonlinear dynamical model; 
the latter, thanks also to the increase of computational facilities, has provided convincing 
evidence of this result by very careful and time consuming numerical simulations. Anyway, 
in both cases strong chaos has been argued to be the basic mechanism responsible for 
normal transport properties. In this perspective, one could conjecture that the doubtful 
results obtained for chaotic potentials, yielding interactions smoother than elastic collisions 
could deserve extremely lengthy numerical simulations before showing standard transport 
properties (e.g., see ||T^Jl3[] ). 

Unfortunately, the explanation of the above mentioned results is quite different. Both 
ding-models can be thought as interacting particles bounded to a local, or substrate, po- 
tential. By substituting the elastic collisions - responsible for the strong chaotic behaviour 
of these models - with any nonlinear nearest-neighbour interaction potential, numerical ev- 
idence of finite thermal conductivity is achieved (for instance, see ref. [|T5). Conversely, 
if the substrate term is absent small wavenumber, i.e. hydrodynainic, modes, as a con- 



sequence of the conservation of momentum , exhibit a peculiar slow relaxation [|15|. This 
effect can be traced back to the vanishing of the generalized frequency uj{k) in the limit of 
small wavenumber, k ^ 0, characterizing the dispersion relation in presence of an acoustic 



band |T^. In practice, the frequency associated with the local potential lifts the acoustic 
band from zero, thus making possible the diffusion of the hydrodynainic modes energy. This 
notwithstanding, it is certainly non-trivial that a purely deterministic chaotic dynamics can 
give rise to an effective diffusive behaviour. Nonetheless, chaos alone is a necessary but not 
sufficient ingredient for obtaining normal trasport properties Q 

Some progress in the understanding of transport properties in nonlinear lattices has 
been made recently by exploiting the analogy with the mode-coupling theory, that applies 
to models of dense fluids rather than diluted gases. Numerical simulations carried out for 
a Id lattice with quartic nonlinearity - the Fermi- Past a- Ulam (FPU) /3-model - have shown 
a power law dependence of the thermal conductivity on the system size: k ~ N°', with 
a = 0.39 ± 0.02 ^^. This result was obtained for a chain in contact with thermal baths 



and by equilibrium measurements of n based on the Green-Kubo formula of linear response 
theory (for technical details , see also Sec. 4). Taking advantage of the results reported in 
|T^, it was possible to conclude that, even for strongly chaotic dynamics, the amplitudes of 
low wavenumber modes evolve like weakly damped stochastic harmonic oscillators, whose 
frequencies are renormalized w.r.t. the standard harmonic component of the FPU /3-model 
by leading nonlinear terms ||18|. On the other hand, large wavenumber modes behave like 
fast variables, governed by a thermal behaviour. This is the typical scheme that applies to 



* We want to point out that normal thermal conductivity has been recently observed also in a class 
of Id hamiltonian models without any substrate potential and intercating via nearest-neighbour 



periodic and double-well chaotic potentials |17 |. 



dense fluids in Fourier space. Self-consistent mode coupling theory (SMT) |T9| , |20[| allows one 
to determine the dependence of k on A^ from the integral of the time-correlation function of 
the total heat flux, yielding in the Id case the power law 

K ~ iV2/5 (3) 

in very good agreement with numerics. This power-law has been recovered also in many 



other models of nonlinearly coupled oscillators like the diatomic Toda chain pT] , p2| , the 
Lennard- Jones (LJ) 6/12 model and the Morse potential [^. All of these models share one 
common feature: they are single-well confining potentials. Also in the light of the results 
recently obtained in Jl^, the scaling law (|^), typical of this class of models, seems to be 
directly related to this property. It is also worth stressing that, due to its generality, SMT 
can be extended to determine the dependence of transport coefficients on the system size in 
higher space dimensions. For the class of models of interest in this paper, it predicts that 
K should diverge logarithmically with A^ in 2d (see Section 4), while in 3d k is normal, i.e. 
independent of A^. 

The aim of this paper is the numerical and analytical study of the thermodynamic limit 
properties of thermal conductivity n in the 2d case, by taking into account two examples of 
single-well confining potentials: the FPU (3 and the LJ-(6/12) models. 

More precisely, in Section 2 we introduce the models and we also discuss some general 
features of the simplified lattice representation of anharmonic 2d solids used in this paper. 
In Section 3 we present the results concerning the dependence of k on A^ obtained by 
two different numerical approaches, based on two different definitions of k. An explicit 
derivation of the analytical prediction of SMT is worked out in Section 4, where we also 
discuss the presence of a super anomalous divergence of transport coefficients in the weakly 
chaotic regime, that typically characterizes the above mentioned models for sufficiently small 
energy [Q (see also [^,^). Conclusions and perspectives are contained in Section 5. 

II. MODELLING HEAT CONDUCTION IN 2D LATTICES 

Most of the studies on the heat conduction problem in anharmonic lattices has been 
devoted to Id systems. Time-consuming numerical simulations were the main limitation for 
large scale analysis of this problem in 2d. In fact, only a few contributions have been worked 
out. For instance in |2^ the dependence of the thermal conductivity k on the temperature 



T is studied in a 2d triangular lattice of unit mass atoms interacting via a LJ 6/12 potential. 
Numerical data are consistent with the expected classical law k ~ T~^, while the dependence 
of n on the system size is not investigated. An interesting contribution in this direction is 



the paper by Jackson and Mistriotis [Q. The authors compare measurements of k in Id 
and 2d Fermi-Pasta-Ulam (FPU) lattices and they conclude that in both cases there is no 
evidence that the transport coefficient is finite in the thermodynamic limit. More recent and 
extended numerical simulations performed for the 2d Toda-lattice [^] have been interpreted 
in favour of the finiteness of k in this case. 

As usual, numerical analysis alone without any piece of a theory can hardly yield con- 
clusive results. On the other hand, as already stressed in [|1^], the dependence of n on the 



system size cannot be adequately described by Peierls' model: at least in the high tem- 
perature (i.e., classical) limit, the perturbative Umklapp processes cannot account for the 
genuine nonlinear effects that characterize such a dependence. 

The theoretical approach proposed in [jl6| and careful numerical analysis of Id systems 



T8| , |29[| , provide a first coherent explanation of the divergence of k{N) in the thermodynamic 
limit. 

In this paper we want to extend this study to 2d systems. For this purpose, let us 
first describe the adopted lattice model of a homogeneous 2d solid. We consider a square 
lattice made of N^ x Ny equal mass (m) atoms. The equilibrium positions of the atoms 
coincide with the lattice sites, labelled by a pair of integer indices (i, j). Without loss of 
generality, the lattice spacing a can be set equal to unit and the origin of the cartesian 
reference frame can be fixed in such a way that 1 < i < A^^^ and 1 < j < Ny. Accordingly, 
the 2d-vector of equilibrium position, ff-, coincides with {i,j). The short-range character 
of interatomic forces in real solids is simplified by assuming that the atoms interact by a 
nearest-neighbour confining potential V{p), that depends on the relative displacement p with 
respect to the equilibrium distance. The natural ordering of the atoms induced by this kind 
of local interaction allows to identify with the same couple of indices {i,j) the corresponding 
atom in any dynamical configuration. Specifically, the model is described by the general 
Hamiltonian 






^) + ^(1^1+1. - %\) + ^(kVi - %l) 



(4) 



where qijit) = Tijit) — ff^, Tijit) is the instantaneous position vector of the (i, j)-atom 
and Pij{t) is the corresponding momentum vector. The time dependence of the canonical 
variables has been omitted in eq.@. Without prejudice of generality we set also m = 1. 

Since we are interested in studying the thermal conductivity k in this 2d model, we have 
also to define the relevant physical observables, namely the temperature T^ and the heat 
flux Jij. Let us start from the local energy density 



«J 



where 

hij = 2 \Pijf + -^[^ili+ij - %\) + ^iWij - ^»-ijl) 

+ ^(|gVi - %-|) + Vi\q,, - gl,_i|)] (6) 

Assuming that local equilibrium holds we obtain immediately a definition of the local tem- 
perature Tij by applying the virial theorem: 



rj.^ , _ , (k) ^J^. _ , (k) tmi^. /7\ 

*^ ~ V/^Jj ^ (fc) / ~ Wy a (fc) / ^^ 



where k = x,y indicates either the x- or the ^/-component of the corresponding vector 
variables. The formal definition of the time average symbol is (•) = limT-_>oo - Jq • dt. 



In numerical experiments the time-average must be performed over a finite time span , in 
general much larger that the local equilibration time. 
For hamiltonian (^) one can use the simple expression 

T„ = ipf) = (pSf ) = ^'V^^^^ (8) 

The heat flux vector is implicitely defined by the continuity equation 

/i^j(t)+V-4(t) = (9) 

By rewriting this equation in Fourier-space and retaining the leading hydrodynamic contri- 
bution (i.e., by applying the large wavelenght limit) one obtains an explicit expression for 
the components of J 



j.(^) = -la 
-'jj 4" 



fxx ( J^) , (^) \ , fyx ( Jy) , (y) 

Jij \Pij ^Pi+ljJ^Jij [Pij ^ Pi+lj 

f^y (J.^) I „W \ I fyy (Jy) , Jy) 
Jij \Pij ^Pij+i)^Jij \Pij ^ Pij+i 



(10) 



where the components of the local vector forces are defined as follows 

cxy _ dV [\qij+i - qij\] _ dV [\qij+i ~ qij\ 



fXy " ■ l\1tJ-f--L lij u r 



dqff '' dcff 

It is worth defining also the space-time average of the heat flux vector 

I J) = &i:k\ (11) 

that will be widely used in the following Sections. 

Finally, since we are interested to investigate thermodynamic limit properties, we should 
impose periodic boundary conditions [jphc)^ that are the standard choice for reducing as 
much as possible undesired boundary effects on bulk properties. 



III. MODELS AND NUMERICAL EXPERIMENTS 

The numerical simulations that we are going to describe in this Section have been per- 
formed for two different single-well potentials: 

Vi(p) = —(? H — p^ , Fermi — Pasta — Ulam [3 — model (12) 

A R R^ 
V2ip) = —777 ^ H 7 , Lennard — Jones 6/12 m^odel (13) 

The reasons for this twofold choice are the following: 



analyzing heat tansport in the 2d version of the widely studied case of potential Vi, 

verifying the theoretical prediction (see next Section) that nearest-neighbour single- 
well nonlinear potentials exhibit the same kind of dependence of k,{N) in the large N 
limit; 

investigating the possible consequences on transport properties of the slowing-down of 
relaxation below the so called Strong Stochasticity Threshold, that was measured for 
potential V2 in pi . 



It is worth stressing that Vi does not contain any natural energy and length scales; all 
numerical simulations with this potential have been performed with /3 = 10~^ and with an 
integration time step At = 10~^, that is two orders of magnitude smaller than the minimum 
harmonic period (xmm = 7r/-\/2). 

At variance with Vi, V2 is characterized by natural length and energy scales, the equi- 
librium distance tq = (2A/By^^ and the well depth W = B'^/AA, respectively. In order to 
make the two models as close as possible we have chosen the parameters A and B in such a 
way that the coefficients of the second and fourth order terms of the Taylor series expansion 
of V2 around its minimum coincide with those of Vi, thus obtaining tq ~ 25, W ~ 8.6 
and Tmin = t^Tq / {<o\/2B) ^ 2.2. In numerical experiments we have adopted a time step 
At = 5 • 10~^, that guarantees a sufficient sampling of the integration algorithm especially 
when strong nonlinearities are explored by the dynamics. 

Numerical measurements of the thermal conductivity k as a function of the lattice size 
A^ have been performed under nonequilibrium and equilibrium conditions: the former are 
inspired by an ideal experiment for the verification of Fourier's law (^, the latter stem 
from linear response theory (see Sec. 4). Both approaches are theoretically equivalent; 
nonetheless, their comparison on a numerical ground is quite instructive. 

Nonetheless, this kind of numerical simulations are quite heavy, so that any trick for 
saving CPU time is worth to be applied. Since we are interested in investigating thermo- 
dynamic limit properties we have to perform measurements keeping the ratio R = Ny/N^ 
constant. On the other hand, there is a priori no reason for choosing R = 1. In fact, we 
have checked that reliable measurements of k, can be obtained in lattices with i? < 1. As an 
example, in Fig.^ we show the dependence of k, on Ny for N^ = 128 and for potential V2, in 
the presence of thermal baths acting on the boundary atoms at temperatures T^ = 1. and 
Tr = 0.5 (details are given in the next subsection). After a sharp increase for small values 
of Ny thermal conductivity k reaches a plateau already for Ny ^ 20. 

This preliminary analysis has been performed for both potentials Vi and V2 in nonequi- 
librium and in equilibrium simulations: we have concluded that i? = 1/2 is a reasonable 
compromise, valid in all of these cases. 



A. Nonequilibrium measurements of thermal conductivity 

A straightforward way of measuring thermal conductivity n amounts to simulate numer- 
ically a true experiment where the atoms at the left and right edges of the 2d lattice are 
coupled with two thermal baths at different temperatures T^ and Tr. This setting restricts 
the application of pbc to the direction orthogonal to the temperature gradient: we impose 



pbc along the y-axis, while the boundary atoms are coupled to a rigid wall through the force 
link with the missing atom in the x-direction. 

Various models of thermal baths, either stochastic or deterministic ones, are at disposal. 
The results of numerical simulations reported in this paper have been obtained by using the 
Nose-Hoover deterministic model 0,0 f\- It has two advantages with respect to stochastic 
algorithms: it can be easily implemented as an ordinary differential equation and it reduces 
the residual thermal impedence effects at the lattice boundaries. 

The equations of motion are 

Qij = Pij 

dV 
Pij = —g^ - (5i,i + Si^NjdjPij 






1 
1 


\\pl^' 

_2Tl 

'\PnJ^ 


^2 


_ 2Tn 



1 



(14) 



where 6 is the Kronecker symbol. In practice, the integration scheme (|n|) has been imple- 
mented by a standard fourth order Runge-Kutta algorithm. Note that each boundary atom 
is coupled through the momentum equation to its thermal bath variable (, that guaran- 
tees local thermal equilibrium at temparature T^ and Tr on the left and right boundaries, 
respectively. 

Starting from random initial conditions after a sufficiently long transient time stationary 
nonequilibrium evolution is eventually approached. According to the heat equation (||), 
a constant thermal gradient should establish through the lattice in the x-direction and 
(J^^^) > 0, while (J'-^-*) = 0. Here, the symbol (•) indicates the time average over stationary 
nonequilibrium states. The time span, over which good convergence of the time-average is 
obtained, increases with N^: for instance, (9(10^) integration steps are sufficient for A^^^ = 16, 
while for N^ = 128 this time grows up to O{10'^). 

The occurrence of stationary nonequilibrium conditions can be checked by verifying the 
equality 

\{J)\ = {JM) = -{j;rT.(M,\PMj\') , M = l,Ar, (15) 

y j=i 

where the r.h.s. is the average heat flux flowing through the boundaries. Despite in numerical 
simulations finite time-averages never yield exactly (J^^-*) = 0, one finds that |(J)| is very 
well approximated by (J*^^-*). 

Some of the stationary temperature profiles obtained for different values of N^ for po- 
tential Vi and V2 are shown in Figures |^ and |^, respectively. The values of the temperature 
Tj have been averaged also in space over all the Ny atoms with abscissa x = i: 



^ We have checked that other models of thermal baths, in particular stochastic ones, yield the 
same kind of results reported in this paper. 

8 



1 "y 
T^ = ^l:T^, (16) 

It is worth stressing that the local temperature Tij is a time-averaged quantity (see eq.®) 
over the stationary nonequilibrium evolution. 

The lattice length has been rescaled to unit in order to verify the overlap of the profiles 
for increasing values of N^'- in both cases we observe a good data collapse for N^ > 32 
indicating that the temperature gradient in the thermodynamic limit vanishes as 

VT, ~ N-' (17) 

In the FPU /3-model the temperatures of the thermal baths, Tl = 20 and Tr = 10, have 
been chosen for obtaining good ergodic properties of the dynamics, i.e. fast relaxation to 
local equilibrium. Boundary effects of thermal impedence induced by the coupling with the 
thermal baths are reduced for increasing values of N^. The temperature gradient looks quite 
close to a constant, but a more careful inspection shows that the profile has a slight curvature. 
The effect is even more evident for the Lennard- Jones model, whose temperature profiles 
converge to an s-shaped curve. In particular, thermal impedence at the boundaries seems 
to persist also in the large N^ limit, despite the smaller values of T^ = 1 and Tr = 0.5, that 
already guarantee ergodicity for this model. Such deviations from Fourier's law indicate that 
nonlinearities influence boundary effects and also the dependence of n on the temperature. 
On the other hand, we are interested in extracting the scaling of n with the system size 
Nx- It can be obtained on the basis of Fourier's law and eq.(|T^ by plotting the quantity 
(J*^^))iVa. ~ K versus N^ (see Figures ^ and |^). For both 2d models we find evidence of a 
logarithmic scaling of k, with A^^.. 

B. Equilibrium measurements of thermal conductivity 

The above described nonequilibrium measurements of the heat transport coefficient k 
have been made keeping T^ and Tr constant for increasing values of A^^,- The data collapse 
of the temperature profiles implies also that in the thermodinamic limit, A^^. —>■ oo, the 
temperature gradient in the bulk of the chain drops to zero. Accordingly, for increasing 
values of A'^,. nonequilibrium measurements are expected to reproduce better and better linear 
response conditions. The Green-Kubo linear response theory |^ provides an alternative 



definition of k with respect to the one contained in Fourier's law. More precisely, a general 
expression of the heat conductivity tensor for a solid contained in a volume V is given by 
the formula 



V 



K 



fii/ 



KbT^ Jo 



{J^^\t)J^''\0))dt (18) 



where Kb is the Boltzmann constant and the time correlation function of the heat flux 
components along the fi and z/ directions is averaged over equilibrium states at temperature 
T. Note that the symbol (•) indicates now the equilibrium ensemble average. 

In our models of homogeneous 2d nonlinear solids the thermal conductivity corresponds 
indifferently to anyone of the diagonal components of the transport coefficient tensor defined 
in ([18|) . In particular we have chosen to measure 

9 



«: = '^.x = #S lim f\j^^\T)J^^\0))dT (19) 

Kfil t-*oo Jo 



ml 

where J^^^ is the x-component of the space-averaged heat flux vector (pH]). As in nonequi- 
hbrium measurements fixing the ratio R = Ny/N^ = 1/2 is a good compromise for reducing 
CPU time, while maintaining reliable asymptotic estimates of n. 

We have performed numerical simulations at constant energy density, e, by eliminating 
the thermal bath variables C's from eq.(|Hp. Periodic boundary conditions have been imposed 
also in the x-direction, so that the 2d lattice becomes a torus. In practice the integration 
scheme has been implemented by a fourth order Maclachlan-Atela symplectic algorithm ||33|| , 
that is more appropriate than a Runge-Kutta one for this kind of microcanonical simulations. 

The expression on the r.h.s of eq.(^) can be evaluated numerically for finite values of t. 
Accordingly, we define a "finite time" thermal conductivity K{t) by removing from eq.(|18|) 
the limit operation t ^ oo. In order to reduce fluctuations, the heat flux time-correlation 
function has been averaged over a sufficiently large set of random initial conditions, ex- 
tracted from the microcanonical probability distribution. For sufficiently large values of t 
both models (0) and (|T3|) exhibit a logarithmic divergence of K{t). As an example, in Fig- 
ure ^ we report the case of the Lennard- Jones (6-12)-model. Despite highly time consuming 
simulations have been performed for averaging over initial conditions, in this example fluc- 
tuations are still persistent for large values of t. This notwithstanding, numerical data are 
compatible with the expected logarithmic divergence. 



According to the argument reported in [it] this divergence in time can be assimilated to 



a divergence with the system size A^^. In fact, numerical estimates of the time-correlation 
function of the local heat flux, Cii^r) = (^oj(t)</jj(0)) show that in single- well nonlinear 
potentials, due to the "rigidity" of low-wavenumber modes (see next Section), the energy 
propagates with the velocity of sound in a nonlinear medium (see also |T^): 



il + a{E))c (20) 

where c is the velocity of sound due to the harmonic component of the nonlinear potential 
and the renormalization factor a{E) > introduces a dependence of c on the total energy, 
indicating its nonlinear character. We have verified that these features are common to the 
potentials (^ ) and (|T3| ). 

The evaluation of K,{t) for finite time t in the Green-Kubo integral can be assimilated to 
an estimate of the divergence with the system size by the formal relation t = N^/c. In this 
sense, equilibrium simulations confirm the results of nonequilibrium ones. 



IV. TRANSPORT IN STRONG AND WEAK CHAOTIC DYNAMICS 

The conjecture that deterministic chaos is an efficient mechanism for ergodic behaviour 
has been extensively investigated by numerical experiments in many degrees of freedom 
Hamiltonian systems, like those introduced in Section 2. One of the main issues is that, for 
sufficiently high energy density e, the time averages of most observables of thermodynamic 
interest rapidly approach the expectation value predicted by equilibrium ensembles. This 
notwithstanding, below a specific value e^ - the so-called Strong Stochasticity Threshold 

10 



(SST) (see p5|j26|] and the references therein contained) - the equihbration time may increase 
dramatically, despite the persistence of deterministic chaos. It is not our aim, here, to discuss 
why such a slow relaxation mechanism typically occurs in these models. We just want to 
stress that their mild nonlinear character make them quite different from mathematically 
standard chaotic models like K- or A-systems. In fact, the results obtained in |T5| show 
that regularities are present even above the SST. More precisely, in a Id FPU chain the 
amplitudes of low-fc Fourier modes are found to evolve like weakly damped and forced 
harmonic oscillators. Dissipation and forcing epitomize the complex nonlinear mechanism 
of energy exchange among the modes and they are found to vanish in the limit fc — > 0. One 
can argue that these can be only an effect induced by global, i.e. hydrodynamic, conservation 
laws of total momentum and energy. In this sense, such a behaviour is expected to be present 
in any space dimension d, although the possibilities of energy exchanges among the modes 
are expected to become more and more efficient for increasing values of d. These remarks 
suggests that, in a nonlinear model of a solid, Fourier modes can be assimilated to a dense 
fluid rather than to a diluted gas, as in Peierls' phonon theory. 

A. Analytical estimate of anomalous thermal conductivity above the SST 

Relying upon the above conjecture, strongly supported by numerical simulations, we 
want to derive an analytical estimate of the thermodynamic limit behaviour of thermal 
conductivity k for 2d single-well anharmonic lattices. In practice, here we explicitely extend 



to the 2d case the method introduced in |jT6[ and described in detail for the Id case in pA 



The amplitudes of the Fourier modes in a 2d square lattice have the standard expression 



E E gn.ne-^l-'^'"^^'^"J (21) 

y 



N^,N„ m=l n=l 



where g^„ is the canonical space coordinate and the 2d wave-vector k has components 
kx = — ^...^ and ky = — ^ . . . ^. We consider an isolated hamiltonian model of 
the type (^; the equations of motion expressed in terms of Fourier amplitudes and their 
canonically conjugated momenta Pj: read 

Pu = -4Qk + E %' (22) 

kj^k' 

where a;^ obeys the dispersion relation 

(23) 



^1 = 4 



sm — - — h sm —-r- 

iVx Ny 



■^For he sake of simplicity we have set ujq = { ^^ ) 



=ro 
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F^g, is the formal expression of the nonhnear interaction force between modes k and k'. 

Numerical experiments show that \ow-k Fourier modes behave like slow dynamical vari- 
able, whose relaxation time scales are much longer than those of high-fc modes, that play the 
role of fast variables. This analogy with hydrodynamic properties of dense fluids suggests 
the apphcation of SMT |9|,|0l. 

Since the hydrodynamic behaviour is ruled by low-A; modes, the first step amounts to 
project the dynamics onto the subspace of slow variables (Q|, P|) by a proper projection 
operator 



V'X = Y, 



-, 2 


-Ql: 


(XpS*) ^ 


i(^l ) 




(PI ) J 



(24) 



According to linear response theory the equations of motion can be casted in the form |p4|| : 



Qk 



PI 



Pk = -^IQI - 1 m - t')Pl{t')dt' + 4^ 



(25) 



V^)Pf^ is the effective force, that is related by the dissipation- fluctuation 



where i?^ = ( 

theorem to the memory kernel as follows: 

r,^(t) oc (4^(t) ■ 4^(0)) 



(26) 



Note that in the projected equations of motion the harmonic frequency u^ is renormalized 
to the energy dependent frequency 



i^i 



where 



a{E) 



^(l + a{E))u^, 



1 J5{H-E)df 
'P J r'^diH - E)df 



1 



(27) 



(28) 



In the last formula we have introduced the microcanonical average by assuming that equi- 
librium properties are practically guaranteed by the fast thermalization of high-fc modes. As 
observed also in [|T5| the projected slow variables can be interpreted as "nonlinear" Fourier 
modes, whose frequency is renormalized by anharmonic contributions. Numerical simula- 
tions support this approach also thanks to the remarkable quantitative agreement with the 
theoretical prediction (^). 

The projected equations of motion have a clear theoretical interpretation but they are 
practically useless for working out analytical calculations. We have to introduce some sim- 
plifying hypotheses. We suppose that the time scale of slow hydrodynamic variables can be 
unambigously separated form the typical microscopic time scales: specifically, this amounts 
to reduce the memory kernel T to a. 6 distribution and the fluctuating force i? to a stochastic 
force. By introducing an explicit complex- variable representation of the mode amplitudes. 



Q,: 



^k + ^B^k 



12 



we rewrite the set of equations (^) in the approximate form of standard stochastic equations: 



^k+iA+^'A 



B^k + IkBu + ^P, 






(29) 



where 7^ can be interpreted as an effective dissipation coefficient, while fluctuations are now 
assumed to be represented by gaussian white noise processes: 



{%it)ffp{t'))cx6^^,,6{t-t') 
It is convenient to introduce the scalar variable 



(30) 



Wr 



^k ■ B^, 



H ■ Bl 



that satisfies the Langevin type equation: 



Wt. 



-ikW^k + C^ 



(31) 



(32) 



where 



Ut)CAt')) « i\M' + \B^?)ku,Kt - 1') 



According to the definition (|Tl]) the average heat flux vector can be thought as the sum 
of a harmonic term and an anharmonic one, Jh and J4, respectively. The harmonic term 
is obtained by considering only the forces given by the quadratic part of the interaction 
potential: simple calculations show that its components can be written in the form: 



J 



(x) 
H 



J: 



(y) 



H 



k 



J2k^,ky Cky^kyW,: 



^X,Ky 'n.y- n.y 



(33) 



where 



^fc. 



irk. 



sm 



a;,j/ 



N, 



x,y 



Cfc. 



Nx,v duJk, 



71 dk. 



x,y 



In close analogy with the renormalization of the effective frequency of slow variables, it can 
be shown that the leading contribution to J a, stemming from the anharmonic part of the 
interaction potential, is proportional to Jh through an energy (and also model) dependent 
factor C{E), i.e. Ja = C{E)Jh- For instance, the expression of C{E) in the FPU /5-model 



(0) is 






^k 



Bt 



(34) 



Summarizing, the components of the average heat flux vector can be expressed in general 
by the proportionality relations 



J^^^ ^Hk^^Cky^K^k 



(35) 
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Assuming the validity of the generahzed equipartiton theorem 

u;l{\A^f)=ujl{\B^f)=U{E) (36) 

where U{E) is a function of the energy E, and using the solution of (|32D, one can obtain 
an analytical estimate of the heat flux time-correlation function present in the Green-Kubo 
integral 



/'27r PTT rjl^ 

{JuityrHm -U{E)Y.cle-'<i' - I del —--kc\koose)e-^'^^' (37) 

k 

where the last expression is obatined by assuming that in the thermodynamic limit the 
summation can be approximated by an integral. 

The explicit dependence of the dissipation coefficient 7 on the modulus of the wave- vector 
k is provided by SMT [19,^|. Specifically, it predicts that the time-correlation functions of 



hydrodynamic modes decay in time as follows: 

G{k,t) ^ e-""^^^' (38) 

where the ui{k) are complex generalized frequencies ruling the oscillatory and relaxation 
behaviour of the hydrodynamic modes. For 2d homogeneous systems they depend just on 
k = \k\ and their explicit expression in the hydrodynamic limit, fc — > 0, is found to be 

uj{k) =:i ick + vk'^hik (39) 

where c is the velocity of sound in the homogeneous medium and the second addendum of 
the r.h.s. is the dissipation term '^{k) = 7le{uj{k)). Substituting into eg. (|37|) and retaining 
only the leading contribution to the function c{kcos6) in the limit A; — > 0, one obtains the 
estimate 

(^iW/i(0)>oc/;|te-'-~i + o(^) (40) 

Accordingly, the heat conductivity k, is predicted to diverge as In t in the t -^ 00 limit in 2d 
lattices. 

It is worth stressing that SMT predicts also that 'j{k) ~ k^^^ in Id systems, yielding a 
diverging heat conductivity k ~ N'^/^ in the thermodynamic limit. The very good agree- 
ment with numerical experiments (see |TB|) strengthen the conjecture that the anomalous 
properties of transport coefficients in nonlinear lattice models stems from the hydrodynamic 
nature of low-k effective modes. Said differently, total energy and momentum conservation 
laws impose on models with single-well potentials a dispersion relation yielding a subdif- 
fusive behaviour of the energy exchange among these modes and, consequently, ill-defined 
transport coefficients. On the other hand, SMT predicts that in 3d the Green-Kubo integral 
is convergent so that k becomes a well defined quantity. 

This dependence on the space dimension and the generality of predictions (at least for 
what concerns the class of single-well nonlinear potentials) make this piece of theory quite 
elegant and physically sound. 
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B. The effect of weak chaos on anomalus transport 

All the numerical simulations presented in Section 3 have been performed in the strong 
chaotic regime of the dynamics. In this case ergodic behaviour occurs quite rapidly on the 
microscopic time scale, with the exception of hydrodynamic modes, whose slow relaxation, 
due to macroscopic conservation laws, rules the observed divergence of transport coefficients. 
In the above theoretical treatment this corresponds to the assumption that time-correlation 
functions can be estimated on the basis of equilibrium ensemble averages. 

On the other hand, we know that below the SST relaxation to equilibrium may slow- 
down dramatically for most of the physically interesting observables. According to the 
results reported in |2^ the heat flux is expected to belong to this class of observables. Then, 
we expect also that in such a dynamical regime the transport mechanism can be significantly 
modified. 

We have checked this conjecture by considering model (|13D , whose SST has been esti- 
mated p3 at a value of the energy density 



essT = 0.3 (41) 

in the units adopted in this paper. We have performed numerical simulations with thermal 
baths at temperature T^ = 0.1 and Tr = 0.05 for different values of N^ verifying the data 
collapse of the temperature profiles. Following the same approach described in Section 3.1 
we have estimated the dependence of the thermal conductivity on N^. We find evidence of 
a power-law divergence 

K ~ N: (42) 

with a = 0.77 (see Fig.^). Note that below the SST the exponent a is a function of the 
energy density e. 

The power law divergence of k below the SST is confirmed by the numerical estimate of 
the Green-Kubo integral. In Fig. |^ we compare dependence of n on time t for e = 0.1 and 
for e = 3.0 , below and above the SST, respectively. The comparison of the two curves in 
a linear versus logarithmic scale shows that the asymptotic behaviour of K{t) is completely 
different in the two chaotic dynamical regimes. In particular, for e = 3.0 we recover a 
logarithmic divergence, while for e = 0.1 we find a power law divergence K{t) ~ t" with 
a = 0.75. 

For what concerns the interpretation of this last result, one has to observe that for a 
sufficiently long time equilibrium conditions will be eventually approached. The asymptotic 
time dependence of the Green-Kubo integral should turn to the logarithmic growth on any 
finite system. The relevance of this superanomalous effect could be established only by 
evaluating the dependence of the relaxation time on the system size. Estimates obtained for 
relatively small size systems seem to indicate that the crossover time between the power-law 
and the logarithmic behaviour of the Green-Kubo integral icreases more than linearly with 
the system size N^. On the other hand, as we have often pointed out along this paper, 
numerics may can provide crucial insight for the understanding of these phenomena, but 
cannot be assumed as a proof of anything. Since this effect might have very interesting 
physical consequences we hope to work out in the near future an analytical approach apt to 
describe the superanomaly of transport coefficients below the SST. 
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V. CONCLUSIONS AND PERSPECTIVES 

The presence of anomalous thermal conductivity in 2d lattices of atoms coupled by 
nonlinear nearest-neighbour single- well potentials has been verified by numerical experiments 
and analytical estimates. According to SMT, the effect of dimensionality is found to make 
the thermodynamic limit power-law divergence of Id models turn to a logarithmic divergence 
of 2d models, thus confirming the soundness of the general theoretical framework worked 
out to interpret such anomalous properties emerging from strongly chaotic dynamics |T6| , |34 . 



It is also worth mentioning that in this dynamical regime 3d models are predicted to exhibit 
normal transport properties, i.e. finite thermal conductivity. 

In this paper we have also pointed out that a super anomalous behaviour ruled by a power- 
law divergence of the thermal conductivity, seems to characterize transport properties of 2d 
models below the SST. Such a behaviour has been observed also in Id systems where the 
crossing of the SST corresponds to an increase of the power from the value 2/5 towards 
the limit value 1, that is expected to be approached for vanishing energy densities (i.e. in 
the harmonic limit). Here we have decided to report just the results for 2d systems, where 
the comparison between the power-law and the logarithmic divergence stresses the effect of 
the extremely slow relaxation mechanism already observed almost half a century ago in the 
seminal numerical experiment by E. Fermi, J. Pasta and S. Ulam [^. 

At present we are not able to conclude just on the basis of numerical simulations if 
superanomalous transport is a finite size effect or an asymptotic property, that could even 
concern 3d systems. Also in this case, the construction of a suitable theoretical approach 
would greatly help in answering to this interesting question. 
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FIG. 1. The thermal conductivity k versus Ny for the Lennard- Jones potential ([l^ with 
N^ = 128, Tl = 1. and Tr = 0.5 . 
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FIG. 2. Temperature profiles of the 2d FPU-/3 model for different values of N^ 
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FIG. 3. Temperature profiles of the 2d Lennard-Jones 6/12 model model for different values of 
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FIG. 4. Dependence of the thermal conductivity k on the system size N^ for the FPU /?-model; 
Tl = 20. and Tr = 10. . Statistical errors have the size of the symbols and the dotted line is the 
best fit. 
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FIG. 5. Dependence of the thermal conductivity k, on the system size N^ for the Lennard- Jones 
(6/12)-model; T^ = 1- and Tr = 0.5. Statistical errors have the size of the symbols and the dotted 
line is the best fit. 
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FIG. 6. The finite-time thermal conductivity K{t) versus time t for the LJ (6/12)-model with 
Nx = 64 and e = 1.5 . The curve has been obtained by averaging over 100 initial conditons. The 
dashed straight line has been drawn just for comparison with numerical data. 
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FIG. 7. The thermal conductivity n versus the system size A^^^ with Tl = 0.1 and Tr = 0.05 for 
the LJ (6/12)-model. The dotted line is a best fit for the power law n ~ N^ yielding a ~ 0.77ib0.03 . 
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FIG. 8. The finite-time thermal conductivity k versus time t for the LJ (6/12)-model with 
iV^ = 64 below (e = 0.1) and above (e = 3.0) the SST. 
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